Landslide Susceptibility Calculation for Tuscany, Italy¶

This lessons can be run on EarthscapeHub. Click the button below:

Run on EarthscapeHub

Introduction¶

Landslide susceptibility is the likelihood of a landslide occurring on the basis of local terrain conditions to estimate “where” landslides are likely to occur. Here we present a modified Data Component Use Case for Landslide Susceptibility Calculation for north Tuscany to explore how slope cohesion and soil saturation impact susceptibility. This Jupyter Notebook will first demonstrate how to use several CSDMS data components to download topography and soil datasets for an area in Tuscany, Italy then calculate the landslide susceptibility given 1) soil cohesion values and 2) soil saturation scenarios. To translate this conceptual excersice to a real-life scenario, we will later focus on a heavy rainfall event in this area on September 10-12th in 2017.

Map Source: www.mapsofworld.com

Learning objectives¶

Key concepts

  • Generate susceptibility maps for the Tuscany region in Italy

  • Relationship between vegetation and soil properties to slope cohesion

  • Explore the influence of root chosesion (as a proxy for vegetation cover) on the lanscape's vulnerability to landslides

  • Relationship between cohesion, saturation, and failure probability

  • Factor of Safety and Susceptibility

  • Application Question: How might different vegetation scenarios (i.e. post wildfire) impact slope stability?

Skills

  • Get a general overview on how data components are implemented to get topographic and environmental data and run a landslide model

  • Become familiar with a simple configuration that can explore how slopes behave as we change factor od safety parameters

  • Make changes to root cohesion values to explore how susceptibility changes under different conditions.

  • Gain hands-on experience with visualizing output in Python

Classroom organization¶

General¶

In this lab we will explore landslide susceptibility under different vegetation covers in the Tuscany region in Italy. To do this we modified a repository created by Gan et al. 2022 that exemplifies how models can communicate with geospatial data using the CSDMS Data components

You will run this notebook on the CSDMS Jupyterhub using the landslides kernel. You can also run it using HYDROSHARE and the CUAHSI JupyterHub (You will need to have access to both resources).

Please note that as of May 14, 2023 the Pymt library is not updated and you won't be able to run the notebook in your own PC (please contact developers for more information).

More detailed instructions for setup and environment creation the are included in the jupyter notebook from the above repository.

Requirements¶

API keys¶

To use the ERA5 and Topography data componentsyou will need to create API key files to download any dataset from around the world.

Follow the instructions in CDS API Key and Open Topography API Key to get your keys.

When you have your API keys, uncomment the 2 lines in the box below to input them. Once you do it one time you can comment them again so you the code doesn't ask for the keys everytime.

In [ ]:
# from utils import install_api_key
# install_api_key()

General Setup¶

  1. Import the necessary packages for the data components to run
  2. Create folders to organize inputs and outputs There are three folders i.e.
    • configuration file folder
    • cache folder
    • results folder
In [1]:
# import packages 
import os
import warnings
import math
from numpy import inf
import numpy as np
import pandas as pd
import xarray
import xesmf as xe
import rioxarray 
import cftime
from datetime import datetime
from tqdm import trange
import matplotlib.pyplot as plt
from matplotlib import colors
import imageio.v2 as imageio
from IPython.display import Video
from landlab import RasterModelGrid, imshow_grid, imshowhs_grid
from pymt.models import Topography, Era5
import matplotlib.pyplot as plt
from landlab.components import BedrockLandslider, PriorityFloodFlowRouter
from landlab.io import read_esri_ascii

from utils import regrid_data, cal_subsurface_flow_depth, cal_safety_factor

warnings.simplefilter(action='ignore', category=FutureWarning)
plt.rcParams.update({'font.size': 14})
In [2]:
# create folders
study_area = 'italy'

config_dir = os.path.join(os.getcwd(), 'config_files_{}'.format(study_area))
results_dir = os.path.join(os.getcwd(), 'results_{}'.format(study_area)) 
cache_dir = os.path.join(os.getcwd(),'cache_{}'.format(study_area))


for folder in [config_dir, results_dir, cache_dir]:
    if not os.path.isdir(folder):
        os.mkdir(folder)
        print(folder)

Grid data acquisition¶

Data components are used to download and visualize datasets of different nature:

  • Topography: to download the Digital Elevation Model (DEM) data with 90m resolution.
  • ERA5: to download the volumetric soil water and the precipitation data of the study region

Additionally:

  • soil depth data is obtained from SoilGrids
  • Topography data along with RasterModelGrid (Landlab) is used to calculate the slope angle for the study area.
DEM¶

The figure below shows the bounding box of the study area. To change the parameter settings go to 'dem_config.yaml' file.

We quickly realized that this entire area is way to big to run this code, so we will focus on in Lucca and Pistoia

Image Source: Rosi, A., Tofani, V., Tanteri, L. et al. The new landslide inventory of Tuscany (Italy) updated with PS-InSAR: geomorphological features and landslide distribution. Landslides 15, 5–19 (2018). https://doi.org/10.1007/s10346-017-0861-4

Run the following cells to:

  • use the configuration file to initialize a data component
  • use the variable and grid related methods of this data component to get:
    • metadata
    • data values
Pulling DEM from Topography data component¶
In [3]:
# initialize Topography data component
dem = Topography()
dem.initialize(os.path.join(config_dir, 'dem_config.yaml'))
In [4]:
# get DEM variable info
var_name = dem.output_var_names[0]
var_unit = dem.var_units(var_name)
var_location = dem.var_location(var_name)
var_type = dem.var_type(var_name)
var_grid = dem.var_grid(var_name)
var_itemsize = dem.var_itemsize(var_name)
var_nbytes = dem.var_nbytes(var_name)
print('variable_name: {} \nvar_unit: {} \nvar_location: {} \nvar_type: {} \nvar_grid: {} \nvar_itemsize: {}' 
            '\nvar_nbytes: {} \n'. format(var_name, var_unit, var_location, var_type, var_grid, var_itemsize, var_nbytes))
variable_name: land_surface__elevation 
var_unit: degrees 
var_location: face 
var_type: int16 
var_grid: 0 
var_itemsize: 2
var_nbytes: 460800 

In [5]:
# get DEM grid info 
dem_grid_ndim = dem.grid_ndim(var_grid) 
dem_grid_type = dem.grid_type(var_grid)
dem_grid_shape = dem.grid_shape(var_grid)
dem_grid_spacing = dem.grid_spacing(var_grid)
dem_grid_origin = dem.grid_origin(var_grid)

print('grid_ndim: {} \ngrid_type: {} \ngrid_shape: {} \ngrid_spacing: {} \ngrid_origin: {}'.format(
    dem_grid_ndim, dem_grid_type, dem_grid_shape, dem_grid_spacing, dem_grid_origin))
grid_ndim: 2 
grid_type: uniform_rectilinear 
grid_shape: [480 480] 
grid_spacing: [ 0.00083333  0.00083333] 
grid_origin: [ 43.85083333  10.5       ]
In [6]:
# get DEM variable data
dem_data = dem.get_value(var_name)
dem_data_2D = dem_data.reshape(dem_grid_shape)

# get X, Y extent for plot
min_y, min_x = dem_grid_origin
max_y = min_y + dem_grid_spacing[0]*(dem_grid_shape[0]-1)
max_x = min_x + dem_grid_spacing[1]*(dem_grid_shape[1]-1)
dy = dem_grid_spacing[0]/2
dx = dem_grid_spacing[1]/2
dem_extent = [min_x - dx, max_x + dx, min_y - dy, max_y + dy]

# plot DEM data
fig, ax = plt.subplots(1,1,figsize=(10,5))
im = ax.imshow(dem_data_2D, extent=dem_extent)
ax.title.set_text('Topography Data')
ax.set_xlabel('longitude [degree_east]')
ax.set_ylabel('latitude [degree_north]')
fig.colorbar(im,label='elevation(m)')
Out[6]:
<matplotlib.colorbar.Colorbar at 0x7ff43cc9abb0>
In [7]:
# Creating RasterModelGrid of elevation data
elev = dem_data_2D.astype('float')
new_grid = RasterModelGrid(elev.shape, 90.0)
new_grid.add_field("node", "topographic__elevation", np.flipud(elev), clobber=True)
new_grid.imshow("topographic__elevation", grid_units=("deg", "deg"),
    colorbar_label="Elevation (m)")
SoilGrids¶

Soil depth data with 250m resolution from SoilGrids system. In this dataset, the maximum soil depth value is 200cm. Grid with values larger than 200cm represents open water area.

In [8]:
# download soil depth data
soil_raster = rioxarray.open_rasterio("https://files.isric.org/soilgrids/former/2017-03-10/data/BDRICM_M_250m_ll.tif")
soil_depth_data = soil_raster.rio.clip_box(
    minx=10.50,
    miny=43.50,
    maxx=11.25,
    maxy=44.25,
)
In [9]:
## plot soil depth data
soil_depth_data.plot(figsize=(6,5),cbar_kwargs={'label': 'depth(m)'})
soil_depth_data.rio.to_raster(os.path.join(cache_dir,'soil_depth.tif'))
plt.title('Soil Depth')
Out[9]:
Text(0.5, 1.0, 'Soil Depth')
Landlab RMG - Slope¶

Topography data and the RasterModelGrid from Landlab to calculate the slope angle

In [10]:
# calculate slope angle using Topography data
model_grid = RasterModelGrid(dem_data_2D.shape,xy_spacing=(90,90))
slope = model_grid.calc_slope_at_node(elevs=dem_data) # slope in radians, 1D array
slope_angle = slope.reshape(dem_data_2D.shape) # reshape as 2D array
In [11]:
# # plot slope angle
fig, ax = plt.subplots(figsize=(10,5))
im=ax.imshow(slope_angle, extent=dem_extent)
cbar = fig.colorbar(im, label= 'radians')
ax.title.set_text('Slope Angle')
ax.set_xlabel('longitude [degree_east]')
ax.set_ylabel('latitude [degree_north]')
Out[11]:
Text(0, 0.5, 'latitude [degree_north]')

Grid Setup¶

All datasets downloaded or calculated above are regrid in the same grid resolution. Again for specific details refer to original notebook

In [12]:
# define destination grid coordinate using Topography data
dem_y = np.flip(np.arange(dem_grid_shape[0])*dem_grid_spacing[0] + dem_grid_origin[0])
dem_x = np.arange(dem_grid_shape[1])*dem_grid_spacing[1] + dem_grid_origin[1]
dest_coor = {'lon': dem_x,
             'lat': dem_y}
In [13]:
# regrid soil depth data
soil_depth_coor = {'lon': soil_depth_data.x.values,
                   'lat': soil_depth_data.y.values}

soil_depth = regrid_data(soil_depth_data.values[0], soil_depth_coor, dest_coor, 'nearest_s2d') 
soil_depth = soil_depth/100 # units conversion as m
In [14]:
# plot regridded soil depth data
fig, ax = plt.subplots(1,1,figsize=(10,5))
im = ax.imshow(soil_depth, extent=dem_extent)
ax.title.set_text('Regridded Soil Depth Data')
ax.set_xlabel('longitude [degree_east]')
ax.set_ylabel('latitude [degree_north]')
cbar = plt.colorbar(im, label='depth(m)')

Landscape Susceptibility¶

Factor of safety

In geological engineering, it is common to take the ratio of the resisting stresses to driving stresses. This ratio is called the factor of safety (FS). When FS is larger than 1, the slope should be stable, while if it is below 1, the driving stress exceeds the resistance and the slope is likely to fail. FS can be calculated with the following function, and cal_safety_factor( ) is implemented based on this function.

$$ FS = \frac{(C_r + C_s)/h_s\rho_sg}{\sin\theta} + \frac{\cos\theta \tan\phi (1-\frac{h_w}{h_s}\rho_w / \rho_s)}{\sin\theta} $$

where,

  • Cr: root cohesion (Pa kg/ms^2)
  • Cs: soil cohesion (Pa kg/ms^2)
  • hs: soil depth (m)
  • hw: subsurface flow depth (m)
  • ρs: soil density (kg/m^3)
  • ρw: water density (kg/m^3)
  • g: gravity acceleration (m/s^2)
  • θ: slope angle
  • φ: soil internal friction angle

Susceptibility

Susceptibility is the inverse of FS. When the susceptibility is larger than 1, it means that the slope of the area is not stable and susceptible to landslide.

$$ \text{Susceptibility} = \frac{1}{FS} $$

Plot Results

There are four subplots created for a list of different susceptibility values for a specified cohesioin value.

Varying Root Cohesion¶

In [15]:
# Function for plotting raster data grid
def plot_ls(new_grid, i, cr):
    # plt.figure(figsize=(6,5))
    imshowhs_grid(
        new_grid,
        "topographic__elevation",
        drape1 = np.sqrt(new_grid.at_node["susceptibility"]),
        plot_type = "Drape1",
        azdeg = 200,
        altdeg = 20,
        alpha = 0.85,
        # var_name = "LS",
        # var_units = r"m",
        grid_units = ("(m)", "(m)"),
        default_fontsize = 10,
        cbar_tick_size = 10,
        cmap = "hot_r",
        limits = [0.5, 1.8],
        ticks_km = False,
        # cbar_width = "100%",
        # cbar_or = "vertical",
        colorbar_label_y = -55,
        # add_label_bbox = True,
        # bbox_to_anchor = [1.03, 0.3, 0.075, 14],
        thres_drape1 = 0.01)
    plt.savefig(results_dir+f'/root_cohesion_plot_{i}.png')
    plt.show()
In [16]:
# Specify minimum and maximum root cohesion parameters for the area
min_Cr = 3000
max_Cr = 12000

# Do a list with 4 equally spaced values within that range 
step_Cr = int(((max_Cr + min_Cr) - min_Cr) / 4)
root_c = [i for i in range(min_Cr, max_Cr + 1, step_Cr)] #list of values for root cohesion repoorted for area
In [17]:
# define mask for non-data area (water bodies, etc)
mask = (slope_angle == 0) & (soil_depth > 2.0)

# we assume constant saturation so that relative wetness is one
subsurface_flow_depth = soil_depth 

# plot susceptibility 
# fig = plt.figure(figsize=(16,10))
# fig.suptitle("Susceptibility of landslide for different root cohesion values")

#change plot grid size if you want to explore more parameters in the above list
nrows, ncols = 2,2

# calculate FS for each root cohesion value
for i,cr in enumerate(root_c):
    safety_factor = cal_safety_factor(slope_angle, subsurface_flow_depth, soil_depth, 
                                  root_cohesion=cr, soil_cohesion=5000, soil_bulk_density=1300,
                                  soil_internal_friction_angle=35)
    
    # calculate susceptibility
    susceptibility = 1.0 / safety_factor
    susceptibility = np.where(~mask, susceptibility, np.nan) #water bodies
    susceptibility = np.where(susceptibility > 0.5, susceptibility, 0)
    susp = new_grid.add_field("node", "susceptibility", np.flipud(susceptibility), clobber=True)
    ax = fig.add_subplot(nrows, ncols, i+1)
    plt.title('Cr = '+str(cr)+" (Pa kg/ms$^{2}$)", loc = 'center', fontdict={'fontsize': 12})
    plot_ls(new_grid, i, cr)
    # im_data = susceptibility
    
###### Comment the 2 above lines and uncomment below to check the differences""   
#     if cr == 5000: 
#         ax = fig.add_subplot(nrows, ncols, 1)
#         im_data = susceptibility
#         reference_sus = np.copy(susceptibility) 
        
#     else: 
#         ax = fig.add_subplot(nrows, ncols, i+1)
#         im_data = susceptibility - reference_sus

    # im_sus = ax.imshow(im_data, cmap='BrBG_r', extent=dem_extent)
    # plt.colorbar(im_sus, ax=ax)
    # im_sus.set_clim(0,1.1)
    # plt.tight_layout()

What differences can you notice in the landscape?

Remember that a susceptibility close to 1 means the landscape is prone to landslides and when it's close to zero it is stable.

Varying root cohesion AND soil saturation¶

Now, we will see how our maps change as we vary two parameters in the Factory of Safety equation, root cohesion and soil saturation

Here we look at three (3) different saturation level scenarios (0%, 50%, and 100%)

Now, what differences can you notice in the landscape with different saturation levels?

Again, remember that a susceptibility close to 1 means the landscape is prone to landslides and when it's close to zero it is stable.

In [18]:
# Specify minimum and maximum root cohesion parameters for the area
min_Cr=3000
max_Cr=12000

# Do a list with 4 equally spaced values within that range 
step_Cr=int(((max_Cr+min_Cr)-min_Cr)/3)
root_c=[i for i in range(min_Cr,max_Cr+1,step_Cr)] #list of values for root cohesion repoorted for area
#root_c=[3000, 8000, 12000] 

#sats=[0.1, 0.5, 1.0]

# define mask for non-data area (water bodies, etc)
mask = (slope_angle==0)&(soil_depth>2.0)

# plot susceptibility 
fig = plt.figure(figsize=(12,10))

################# Want to put x axis label at top, want tick marks centered on the columns #####################
################# Want y axis ticks centered on rows ###########################################################

# fig.suptitle("Susceptibility of landslide", y=0.95) #for different soil saturations and root cohesion values")

#change plot grid size if you want to explore more parameters in the above list
nrows, ncols = 3,3 #rows are cohesion values, columns are saturation

gs = fig.add_gridspec(nrows=nrows+2, ncols=ncols+1, width_ratios=[0.1, 1, 1, 1], height_ratios=[0.1, 1, 1, 1, 0.1], wspace=0.2,hspace=0.3)
xlab_ax = fig.add_subplot(gs[0,:])
ylab_ax = fig.add_subplot(gs[1:-1,0])

xlab_ax.text(x=0.5, y=0.05, ha='center', s='Saturation %')
ylab_ax.text(x=0.05, y=0.5, rotation=90, va='center', s='Root Cohesion (Pa kg/ms^2)')

xlab_ax.axis('off')
ylab_ax.axis('off')

w=0

percent = np.arange(0,11,5)
sat_depth = [p/10 * soil_depth for p in percent]

for i, cr in enumerate(root_c, start=1):
    
    for j, subsurface_flow_depth in enumerate(sat_depth, start=1): # Only interested in 0%, 50%, and 100% Saturation
        
        # if j==0:
        #     subsurface_flow_depth = 0
            
        safety_factor = cal_safety_factor(slope_angle, subsurface_flow_depth, soil_depth, 
                                      root_cohesion=cr, soil_cohesion=5000, soil_bulk_density=1300,
                                      soil_internal_friction_angle=35)

        # calculate susceptibility
        susceptibility = 1.0 / safety_factor
        susceptibility = np.where(~mask, susceptibility, np.nan) #water bodies
        
        # Changing the index 
        w=w+1
    
        # plot 
        # ax = fig.add_subplot(nrows,ncols,w)
        ax = fig.add_subplot(gs[i,j], box_aspect=1)
        im_data = susceptibility
        
        ax.tick_params(axis='both', which='major', labelsize=8)
        ax.tick_params(axis='both', which='minor', labelsize=8)
        im_sus = ax.imshow(im_data, cmap='BrBG_r', extent=dem_extent)
        
        if j == 1:
            ax.set_ylabel(cr)
        if i == 1:
            ax.set_xlabel(10*percent[j-1])
            ax.xaxis.set_label_position('top')

        # plt.colorbar(im_sus, ax=ax)
        im_sus.set_clim(0,1.0)

        # plt.close(fig)
cb_ax = fig.add_subplot(gs[-1, :])
cbar = fig.colorbar(im_sus, cax=cb_ax, location='bottom', orientation='horizontal')
cbar.set_label('Susceptibility of Landslide')
# im_sus.set_clim(0,1.1)
# fig.savefig(os.path.join(results_dir, 'sus_root_{}.png'.format(j)))

Let's overlay this on a Hillshade of the area to better visualize¶

This process is similiar to what we did for the original root cohesion comparison

In [19]:
# Function for plotting raster data grid
def plot_ls(new_grid, i, j):
    # plt.figure(figsize=(6,5))
    imshowhs_grid(
        new_grid,
        "topographic__elevation",
        drape1 = np.sqrt(new_grid.at_node["susceptibility"]),
        plot_type = "Drape1",
        azdeg = 200,
        altdeg = 20,
        alpha = 0.85,
        # var_name = "LS",
        # var_units = r"m",
        grid_units = ("(m)", "(m)"),
        default_fontsize = 10,
        cbar_tick_size = 10,
        cmap = "hot_r",
        limits = [0.5, 1.8],
        ticks_km = False,
        # cbar_width = "100%",
        # cbar_or = "vertical",
        colorbar_label_y = -55,
        # add_label_bbox = True,
        # bbox_to_anchor = [1.03, 0.3, 0.075, 14],
        thres_drape1 = 0.01)
    plt.savefig(results_dir+f'/sus_root_plot_{i}_{j}.png')
    plt.show()
In [20]:
# Specify minimum and maximum root cohesion parameters for the area
min_Cr=3000
max_Cr=12000

# Do a list with 4 equally spaced values within that range 
step_Cr=int(((max_Cr+min_Cr)-min_Cr)/3)
root_c=[i for i in range(min_Cr,max_Cr+1,step_Cr)] #list of values for root cohesion repoorted for area
#root_c=[3000, 8000, 12000] 

#sats=[0.1, 0.5, 1.0]

# define mask for non-data area (water bodies, etc)
mask = (slope_angle==0)&(soil_depth>2.0)

# plot susceptibility 
# fig = plt.figure(figsize=(16,16))
# plt.xlabel('Saturation %')
# plt.ylabel('Root Cohesion (Pa kg/ms^2)')
# cb_ax = fig.add_axes([1.02, 0.1, 0.02, 0.8])
# cbar = fig.colorbar(im_sus, cax=cb_ax)
# im_sus.set_clim(0,1.1)

################# Want to put x axis label at top, want tick marks centered on the columns #####################
################# Want y axis ticks centered on rows ###########################################################

fig.suptitle("Susceptibility of landslide") #for different soil saturations and root cohesion values")

#change plot grid size if you want to explore more parameters in the above list
nrows, ncols = 3,3 #rows are cohesion values, columns are saturation
    
w=0

for i,cr in enumerate(root_c):  
    
    for j in range(0,11,5): # Only interested in 0%, 50%, and 100% Saturation
        
        if j==0:
            subsurface_flow_depth = 0
            
        subsurface_flow_depth = soil_depth*j/10
        safety_factor = cal_safety_factor(slope_angle, subsurface_flow_depth, soil_depth, 
                                      root_cohesion=cr, soil_cohesion=5000, soil_bulk_density=1300,
                                      soil_internal_friction_angle=35)

        # calculate susceptibility
        susceptibility = 1.0 / safety_factor
        susceptibility = np.where(~mask, susceptibility, np.nan) #water bodies
        
        susceptibility = np.where(susceptibility > 0.5, susceptibility, 0)
        susp = new_grid.add_field("node", "susceptibility", np.flipud(susceptibility), clobber=True)
        # ax = fig.add_subplot(nrows, ncols, i+1)
        plt.title('Cr = '+str(cr)+" (Pa kg/ms$^{2}$)", loc = 'center', fontdict={'fontsize': 12})
        plot_ls(new_grid, i, j)

ERA5¶

If the model is run in its complete form:

  • volumetric soil water data will be used for calculating the susceptibility

  • precipitation data will help visualize

The 'era5_config.yaml' file includes the parameter settings of this data component.

Run the following cells to:

  • initialize ERA5
  • see how to specify grid and time methods to get metadata and data values
In [21]:
# initialize ERA5 data component
era5 = Era5()
era5.initialize(os.path.join(config_dir,'era5_config.yaml'))
2023-05-17 04:42:09,421 INFO Welcome to the CDS
2023-05-17 04:42:09,421 INFO Sending request to https://cds.climate.copernicus.eu/api/v2/resources/reanalysis-era5-single-levels
2023-05-17 04:42:09,809 INFO Request is queued
2023-05-17 04:42:10,982 INFO Request is running
2023-05-17 04:42:31,644 INFO Request is completed
2023-05-17 04:42:31,645 INFO Downloading https://download-0007-clone.copernicus-climate.eu/cache-compute-0007/cache/data8/adaptor.mars.internal-1684298544.1200817-10618-3-07c501a3-bafd-41b3-87de-3675db8f0d75.nc to ./cache_italy/era5.nc (6.2K)
2023-05-17 04:42:32,488 INFO Download rate 7.4K/s
In [22]:
# get ERA5 variable info
for var_name in era5.output_var_names:
    var_unit = era5.var_units(var_name)
    var_location = era5.var_location(var_name)
    var_type = era5.var_type(var_name)
    var_grid = era5.var_grid(var_name)
    var_itemsize = era5.var_itemsize(var_name)
    var_nbytes = era5.var_nbytes(var_name)
    # print('variable_name: {} \nvar_unit: {} \nvar_location: {} \nvar_type: {} \nvar_grid: {} \nvar_itemsize: {}' 
    #         '\nvar_nbytes: {} \n'. format(var_name, var_unit, var_location, var_type, var_grid, var_itemsize, var_nbytes))
In [23]:
# get ERA5 grid info
era5_grid_ndim = era5.grid_ndim(var_grid) 
era5_grid_type = era5.grid_type(var_grid)
era5_grid_shape = era5.grid_shape(var_grid)
era5_grid_spacing = era5.grid_spacing(var_grid)
era5_grid_origin = era5.grid_origin(var_grid)

# print('grid_ndim: {} \ngrid_type: {} \ngrid_shape: {} \ngrid_spacing: {} \ngrid_origin: {}'.format(
#     era5_grid_ndim, era5_grid_type, era5_grid_shape, era5_grid_spacing, era5_grid_origin))
In [24]:
# get ERA5 time info
era5_start_time = era5.start_time
era5_end_time = era5.end_time
era5_time_step = era5.time_step
era5_time_unit = era5.time_units
era5_time_steps = int((era5_end_time - era5_start_time)/era5_time_step) + 1

# print('start_time:{} \nend_time:{} \ntime_step:{} \ntime_unit:{} \ntime_steps:{}'.format(
#     era5_start_time, era5_end_time, era5_time_step, era5_time_unit, era5_time_steps))
In [25]:
# get ERA5 variables data and plot (at the first time step)
fig = plt.figure(figsize=(18,16))
nrows, ncols = 3, 2
i = 1

for var_name in era5.output_var_names:
    ax = fig.add_subplot(nrows, ncols, i)
    var_unit = era5.var_units(var_name)
    
    # get variable data    
    era5_data = era5.get_value(var_name)
    era5_data_2D = era5_data.reshape(era5_grid_shape)
    
    # get X, Y extent for plot
    min_y, min_x = era5_grid_origin
    max_y = min_y + era5_grid_spacing[0]*(era5_grid_shape[0]-1)
    max_x = min_x + era5_grid_spacing[1]*(era5_grid_shape[1]-1)
    dy = era5_grid_spacing[0]/2
    dx = era5_grid_spacing[1]/2
    era5_extent = [min_x - dx, max_x + dx, min_y - dy, max_y + dy]

    # plot data
    im = ax.imshow(era5_data_2D, extent=era5_extent, cmap='Blues')
    ax.title.set_text('{} ({})'.format(var_name,var_unit ))
    ax.set_xlabel('longitude [degree_east]')
    ax.set_ylabel('latitude [degree_north]')
    cbar = plt.colorbar(im, ax=ax)
    plt.tight_layout()
    i += 1

References¶

  • Calvello, M., Pecoraro, G. FraneItalia, (2018), A catalog of recent Italian landslides. Geoenviron Disasters 5, 13. https://doi.org/10.1186/s40677-018-0105-5
  • Gan, T., Campforts, B., Tucker, G. E., Overeem, I. (2022). Data Component Use Case for Landslide Susceptibility Calculation, HydroShare, https://www.hydroshare.org/resource/df5fa2f5d1b74be4bf0a049e1e59889c/
  • Montgomery, D. R., and Dietrich, W. E. (1994), A physically based model for the topographic control on shallow landsliding, Water Resour. Res., 30( 4), 1153– 1171, https://doi.org/10.1029/93WR02979.
  • Rosi, A., Tofani, V., Tanteri, L. et al.,(2018), The new landslide inventory of Tuscany (Italy) updated with PS-InSAR: geomorphological features and landslide distribution. Landslides 15, 5–19 (2018), https://doi.org/10.1007/s10346-017-0861-4.
  • Schwarz, M., Preti, F., Giadrossich, F., Lehmann, P., Or, D. (2010). Quantifying the role of vegetation in slope stability: A case study in Tuscany (Italy), Ecological Engineering, Volume 36, Issue 3, Pages 285-291, ISSN 0925-8574,https://doi.org/10.1016/j.ecoleng.2009.06.014.
  • Strauch, R., Istanbulluoglu, E., Nudurupati, S. S., Bandaragoda, C., Gasparini, N. M., and Tucker, G. E. (2018), A hydroclimatological approach to predicting regional landslide probability using Landlab, Earth Surf. Dynam., 6, 49–75, https://doi.org/10.5194/esurf-6-49-2018